import sys
sys.path.append("../reading_data")

import logging
from logging import debug
logging.basicConfig(format='%(levelname)s:%(message)s', level=logging.DEBUG)

import datetime

import numpy as np
from scipy import interpolate

import matplotlib.pyplot as plt
from matplotlib.dates import num2date
from matplotlib.dates import date2num


from read_from_h5 import load_data


def get_cp(beta, rate):
    Cp_tab = np.genfromtxt("CpTbl.txt")
    Tp_speed_rate = np.genfromtxt("TpSpdRatio.txt")
    Ptch_angl = np.genfromtxt("PtchAngl.txt")
    f = interpolate.RectBivariateSpline(Tp_speed_rate,
                                        Ptch_angl,
                                        Cp_tab,
                                        kx=1,
                                        ky=1,
                                        s=0)
    z = np.array([f(rat, bet)[0, 0] for rat, bet in zip(rate, beta)])
    return z

def get_cp_from_func(beta, rate):
    z = -3.2160923889e-001 + \
        -8.7961435621e-003 * rate ** 2 + \
        -8.8126938607e-004 * beta ** 2 + \
        -2.1520547106e-003 * rate * beta + \
        3.9836543835e-004 * rate * beta ** 2 + \
        1.6132999265e-004 * rate ** 2 * beta + \
        -8.1022625082e-005 * rate ** 2 * beta ** 2 + \
        5.7359338369e-008 * rate ** 3 * beta ** 3 + \
        1.6304486750e-001 * rate + \
        1.0967507910e-002 * beta + \
        -2.1520547106e-003 * rate * beta
    return z

def plot_pa_pe():
    """
    >>> plot_pa_pe()
    """
    #print [get_cp_from_func(0, i) for i in range(1,20)]
    turbid = "WH1108"
    h5_file = "../../h5data/wh1.h5"
    filt_strings = ["(activepoweravg > 0)",
                    "(activepoweravg < 2200)",
                    "(time >= " + str(date2num(datetime.datetime(2000, 01, 01))) + " )",
                    "(time <= " + str(date2num(datetime.datetime(2012,01,01))) + " )"]
                    #"(time <= " + str(date2num(datetime.datetime(2012, 01, 01))) + " )",
                    #"(turbid == '" + turbid + "' )"]
    debug(filt_strings)

    data = load_data(h5_file, filt=filt_strings)
    rot_vec = data["rotorrpmavg"]
    pitch_vec = data["bladespitchangleavg"]
    wind_vec = data["windspeedavg"]
    time_vec = data["time"]
    power_vec = data["activepoweravg"]
    winddir_vec = data["winddirabsoluteavg"]
    cosphi_vec = data["cosphiavg"]
    power_vec *= 1000
    rho = 1.24
    R = 40
    pi = np.pi
    rpm_to_rad = pi / 30
    rec_fac = 0.9
    mean_wind = 10
    v_old = wind_vec
    v = wind_vec
    v = wind_vec + (1 - rec_fac) * (wind_vec - mean_wind)
    #v_old = v
    w_r = rot_vec

    rate = R * (rpm_to_rad * w_r) / v

    beta = pitch_vec
    #beta = 2 * np.ones(len(pitch_vec))
    Cp = get_cp(beta, rate)
    Pa = 0.5 * rho * v ** 3.0 * pi * (R ** 2.) * Cp
    #Pa += 100000
    Pa[Pa < 0] = 0
    fig, ax = plt.subplots(1, 1)
    #ax.plot(num2date(time_vec), pow_diff, 'o', alpha=0.05)
    #WINDOW = 100
    #weightings = np.repeat(1.0, WINDOW) / WINDOW
    #wieght = np.convolve(pow_diff, weightings)[:-(WINDOW-1)]
    #ax.plot(time_vec, wieght, 'o')
    #ax.plot(wieght)
    #ax.plot(num2date(time_vec), Pa, 'o', label="Pa")
    #ax.plot(num2date(time_vec), power_vec, 'o', label="Pe")
    #ax.plot(v_old, Pa, 'bo', alpha=0.004)
    ax.plot(v_old, power_vec, 'go', alpha=0.004, markersize=3)
    #ax.plot(min(v_old), min(Pa), 'bo', label="Pa", alpha=1)
    ax.plot(min(v_old), min(power_vec), 'go', label="Pe", alpha=1)
    #ax.plot(pitch_vec, Pa, 'o', label="Pa", alpha=0.01)
    #ax.plot(pitch_vec, power_vec, 'o', label="Pe", alpha=0.01)
    #ax.plot(wind_vec, pitch_vec, 'ro', label="pitch", alpha=0.01)
    #ax.plot(wind_vec, rot_vec, 'bo', label="rotor", alpha=0.01)
    plt.xlim(0, 25)
    plt.ylim(0, 2010000)
    #ax.legend()
    #fig, ax = plt.subplots(1,1)
    #ax.plot(num2date(time_vec), pow_diff, 'o', label="diff")
    #ax.plot(pitch_vec, Pa, 'o', label="Pa")
    #ax.plot(pitch_vec, power_vec*1000, 'o', label="Pe")
    #ax.plot(wind_vec, Cp, 'o', label="Pe", alpha=0.05)
    ax.legend()
    fig.autofmt_xdate()
    plt.show()


if __name__ == '__main__':
    import energi_balance
    import doctest
    doctest.testmod(energi_balance, verbose=False)
